Localization by target motion analysis (TMA) is a standard technique used in passive sonar systems. There are many TMA algorithms available at present which use as their inputs a time sequence of bearing measurements and, if the target is narrowband, frequency measurements as well. A discrete sequence of parameter estimates must be obtained at the output of the sonar processor before performing the target motion analysis. This parameter estimation, in its simplest form, reduces to measuring the coordinates of a peak in a two-dimensional frequency-azimuth (FRAZ) power spectrum. It is desirable that the TMA process be automated with the data extraction being performed without operator intervention. In this manner, the bearing/frequency history of multiple targets can be obtained without dedicating an operator to the task of data extraction. However, the automatic following of a spectral peak as time progresses poses difficulties, particularly for passive sonars.
There are several reasons why target motion analysis with passive sonar is inherently difficult. The primary difficulty is the variation of the received acoustic signal resulting from the complex propagation structure of the oceanic medium through which it transits. The received signal may, for instance, have transited through a spectral region of elevated broadband noise produced by a noisy surface vessel. It is, as a result, quite common for signals to fade in and out and at least some fades may last for several integration periods of the sonar processor. The effect of these fades, for whatever cause, is a loss of signal-to-noise ratio (SNR) that can periodically place the signal below the detection threshold.
This sort of signal fading is an impediment to the successful implementation of any automatic signal-follower that attempts to follow a signal by associating a peak in each incoming spectrum with a peak in the previous spectrum. This requires that some method be found to "ride out" fades and associate a reappearing spectral peak with one that existed several integration periods previously. Furthermore, as the SNR of a received signal fluctuates, a noise peak in the immediate vicinity of the signal peak may exceed it in magnitude yielding spurious values for the estimated bearing and frequency. Therefore, it is necessary to pre-filter the raw parameter estimates before they are used as inputs to the TMA algorithm in order to eliminate outliers, i.e. values so far removed from others that their presence may adversely affect the TMA algorithm. However, this pre-filtering stage is usually not robust so the results cannot be confidently used for target motion analysis without examining them and, if necessary, manually correcting them.
A method of target motion analysis according to the present invention can deal with many of the above-mentioned difficulties. This method can be used for both bearings-only and Doppler/bearing TMA. The method, according to the present invention, integrates raw FRAZ spectra in both time and space along the coordinate trajectory generated by a hypothesized target track. That is, frequency and bearing values (the coordinates in frequency-azimuth space) that would be observed by the sonar system are computed for each point on a hypothesized target track corresponding to an integration period of the sonar processor. A long-term spectral value for the entire observation time is then obtained by non-coherently time integrating the spectral values along the coordinate trajectory generated by the hypothesized track. This is in contrast to known sonar signal processing in which the time integration is always for the same frequency bin and the same beam or azimuth.
The search for a target is conducted by integrating over a large number of hypothesized constant-velocity tracks. The search is organized by constructing a dense grid of start and end positions connected by straight-line tracks. The speed required for the target to traverse each hypothesized track is readily found since the total observation time is known. In the narrowband case, when Doppler effects are present, the rest frequency of the emitted tonal must be known in order to compute the frequency that would be observed by the sonar as a function of time. However, the rest frequency is generally unknown and this makes it necessary for the search to proceed over a discrete set of hypothesized rest frequencies as well. The hypothesized track and rest frequency that produce the maximum integrated spectral value are taken as the estimated parameters for the target at the end of the search.
The above-described method avoids many of the difficulties which occur in the implementation of TMA algorithms. There is no need to extract discrete estimates of the target bearing and frequency at the sonar output since raw spectral data is used as the inputs. Furthermore, the algorithm always integrates along a fixed coordinate trajectory corresponding to a hypothesized target track so that a low SNR or fading target is correctly integrated over the total observation time. The algorithm does not "lock on" to a strong noise peak or "lose lock" during a signal fade whereas a recursive signal follower might "lock on" or "lose lock" under those conditions. The method not only provides an estimate of the target range and course but also, as a consequence of its long-term integration of raw spectral estimates, provides an increased detection capability as well. It should be noted that the total observation time may be on the order of hours whereas the integration time of a typical sonar processor is in the range of 5-80 seconds for each spectral output of the sonar processor, i.e. for each line on a spectrogram display. The overall non-coherent integration time possible in a conventional sonar processor is set by the size of the display surface that the sonar operators may view and the time that the target tonal remains in a single beam and frequency bin. The latter constraint on integration time is the determining factor for the narrow beams and frequency bins produced by a modern towed-array signal processor. This constraint is removed by the TMA method according to the present invention since the integration varies in both the bearing and frequency coordinates.
Target motion analysis (TMA) can be performed directly from the beam spectra produced by a towed line array sonar processor using an algorithm developed at the Defence Research Establishment Pacific. The data available to the algorithm is taken from the output of the front-end sonar processor, consisting of a beamformer followed by a spectrum analyzer. In the algorithm, B.sub.k (f,cos.beta.) denotes the power spectrum at a frequency f and beam angle .beta. at time index k, where 1.ltoreq.k.ltoreq.K. The dependence of B.sub.k on the steering angle .beta. is expressed through cos.beta. and initially f and cos.beta. can be allowed to assume arbitrary (not necessarily discrete) values for convenience. B.sub.k can be thought of as a two-dimensional frequency-azimuth (FRAZ) spectrum, the spectrum B.sub.k being referred to simply as a FRAZ. The integration time used in producing each FRAZ at a particular time k is short, i.e. typically 15-120 seconds which is the "update time" of the front-end processor. The idea of this TMA method is to form long-term integrated spectral values from the short-term values according to the equation ##EQU1## where f.sub.k and B.sub.k are the frequency and angle coordinates that would theoretically be observed by the sonar at time k from an assumed target. The long-term integrated power in Equation (1) should attain a maximum when the coordinate trajectory (f.sub.k,cos.beta..sub.k) of an assumed target aligns with the coordinate trajectory of an actual target. It then becomes apparent that TMA based on an exhaustive search is feasible by searching over a large number of hypothesized coordinate sequences (f.sub.k,cos.beta..sub.k) and looking for peaks in the long-term integrated power.
The search over a large number of hypothesized coordinate sequences presents a computational problem in calculating the coordinate trajectories (f.sub.k,cos.beta..sub.k) and performing the addition in Equation (1) along the trajectories. This computational burden is sufficiently great that specialized hardware becomes a necessity when the number of trajectories is large. This also requires that a systematic method of generating the hypothesized coordinate trajectories be chosen so as to represent the majority of potential targets without undue computational complexity.
In order to minimize the computations needed to generate a target's coordinate trajectory only constant velocity tracks are searched since this will lessen the number of parameters required to describe a target. The coordinate vector of a target relative to the origin of a coordinate system at integration time k is denoted by r.sub.k =(X.sub.k,Y.sub.k)so that r.sub.1 is the target's initial position and r.sub.K its final position. A Doppler/bearing TMA is possible if the target is emitting a narrowband tonal, in which case f.sub.o would denote the rest frequency of the tonal. The vector of parameters used to specify a target is denoted by p so that the parameter vector p for bearings-only TMA will be p=(r.sub.1,r.sub.K) whereas for Doppler/bearing TMA, it will be p=(r.sub.1,r.sub.K,f.sub.o). Since the bearings-only case is computationally simpler than the Doppler/bearing case, the latter will be considered in the following description with comments being inserted where needed to illustrate the differences. Assuming constant velocity tracks, the coordinate trajectory (f.sub.k,cos.beta..sub.k) is determined completely from the vector p. The long-term integrated spectral power is denoted by &lt;B(p)&gt; so that Equation (1) can be written as: ##EQU2## where f.sub.k =f.sub.k (p) and .beta..sub.k =.beta..sub.k (p). A search can then be conducted over values of p within a pre-defined volume and the value p=argmax &lt;B(p)&gt; that maximizes &lt;B(p)&gt; is taken to be an estimate of the target parameters. The maximizing value p is that value for which (f.sub.k,cos.beta..sub.k) is aligned with a spectral peak most consistently during the long-term integration. It should be noted that strong peaks are automatically weighted more heavily than weaker peaks. If the search volume contains multiple targets, then local maxima as well as a global maximum must be found.
The coordinate trajectory (f.sub.k,cos.beta..sub.k) can now be computed given the parameter vector p. First of all, it is necessary to distinguish between the beam angle .beta. and the true target bearing .theta. measured with respect to the array heading when the dominant arrival path is a bottom-bounce path. These two angles are connected by the equation cos.beta.=cos.theta.cos.psi. for a towed line array, .psi. being the vertical arrival angle at the array. The depths of the sonar array and target are usually small compared to the water depth H. Therefore, to a good approximation, cos.psi. is given by ##EQU3## where R is in the range to the target and l=0, 1, 2. . . is the arrival order (l=0 for the direct path l=1 for the first bottom-bounce arrival, etc.). The beam angle is equal to the true target bearing for a direct-path arrival when this approximation is used.
Information is also required concerning the track of the towed array in addition to the vector p of target parameters. A two-dimensional vector giving the x,y coordinates of the towed array with respect to the origin at time k is denoted by a.sub.k. Generally, the origin of the coordinate system will be near the array track such as at a.sub.1, the initial array position. The velocity vector of the array at time k is denoted by a.sub.k. Both a.sub.k and a.sub.k are measured quantities which are known before the search begins. An array heading vector of unit magnitude is defined, for convenience, by h.sub.k =a.sub.k /.sub..parallel. a.sub.k.parallel.. The velocity v of the hypothesized target will then be v=(r.sub.K -r.sub.1)/T where T denotes the total observation time. Therefore, the position of the target relative to the array at time k is EQU R.sub.k =r.sub.1 +t.sub.k v-a.sub.k ( 4)
where t.sub.k is the time of the k.sup.th update. Letting R.sub.k =.sub..parallel. R.sub.k.parallel. and defining a unit vector u.sub.k =R.sub.k /R.sub.k that points from the array to the target results in having EQU cos.theta..sub.k =u.sub.k .multidot.h.sub.k ( 5)
and ##EQU4## The angle coordinate cos.sub..beta.k observed at the time k can be computed from Equations (5) and (6) since cos.sub..beta.k =cos.sub..theta.k cos.sub..psi.k. The relative velocity between the target and the array is denoted by v.sub.k .ident.v-a.sub.k and the frequency coordinate at time k is then given by ##EQU5## where C represents the speed of sound in water and is assumed to be constant.
The above-described algorithm requires a search over a volume and the integration specified in Equation (2) must be carried out along the coordinate trajectory (f.sub.k,cos.beta..sub.k) specified by the above equations for each value of p within the volume. It is desirable to frame the algorithm in such a way that it uses only the data stored in the database of a typical sonar processor, i.e. the two-dimensional spectra B(f,cos.beta.) given on a discrete frequency-azimuth grid. Assuming that the sonar processor produces P frequency bins and Q beams with the frequency spacing being constant, as from an FFT, then the P frequencies are given by the equation EQU f.sup.(p) =f.sub.s +p.DELTA.f, 0.ltoreq.p.ltoreq.P-1 (8)
where f.sub.s is the start frequency and .DELTA.f is the spacing between frequency bins. Furthermore, assuming that the Q beams are spaced equally in cos.beta., then EQU cos.beta..sup.(q) =cos.beta..sub.s +q.DELTA.(cos.beta.), 0.ltoreq.q.ltoreq.Q-1 (9)
where cos.beta..sub.S is the cosine of the start angle and .DELTA.(cos.beta.) is the beam spacing. Then given a coordinate trajectory (f.sub.k,cos.beta..sub.k), a sequence of corresponding integer frequency and angle coordinates (p.sub.k,q.sub.k) can be computed such that the differences EQU .vertline.f.sup.(pk) -f.sub.k .vertline. and .vertline.cos.beta..sup.(qk) -cos.beta..sub.k .vertline. (10)
are minimized, and B(f.sub.k,cos.beta..sub.k) is approximated by B(f.sup.(pk),cos.beta..sup.(qk)). In other words, B(f.sub.k,cos.beta..sub.k) is approximated by its value at the closest bin in the pre-computed two-dimensional FRAZ spectrum. This assumes that the FRAZ spectra are computed on a grid that is sufficiently dense that only a small error is incurred by this approximation. A more general approach would use interpolation of several bins from the FRAZ spectrum to get a refined estimate of the spectrum level at the coordinate value (f.sub.k,cos.beta..sub.k). For notational convenience, a matrix B(p,q) is defined by ##EQU6## The long-term integrated sum in Equation (2) is then approximated by ##EQU7##